GL-TR-89-0346 


( 


REMOTE  SOUNDING  OF  ATMOSPHERIC  TEMPERATURE 
PROFILES  USING  THE  OPTICAL  MEASURE  METHOD 


S.  C.  Ou 
K.  N.  Liou 


AD-A244  828 


Center  for  Atmospheric  and  Remote  Sounding  Studies  (CARSS) 

Department  of  Meteorology 

University  of  Utah 

Salt  Lake  City,  Utah  84112 


31  December  1989 


Final  Report 

12  September  1988  -  31  December  1989 


Approved  for  public  release;  distribution  unlimited 


GEOPHYSICS  LABORATORY 

AIR  FORCE  SYSTEMS  COMMAND 

UNITED  STATES  AIR  FORCE 

HANSCOM  AFB,  MASSACHUSETTS  01731-5000 


•8  H  14*008 


92-01226 


This  technical  report  has  been  reviewed  and  is  approved  for  publication. 


Contract  Manager 
Atmospheric  Sciences  Division 


KT  A.  McOATCHEY,  Director 
Cboospheric  Sciences  Division 


This  report  has  been  reviewed  by  the  ESD  Public  Affairs  Office  (PA)  and  is 
releasable  to  the  National  Technical  Information  Service  (NTIS). 


Qualified  requestors  may  obtain  additional  copies  froa  the  Defense  Technical 
Information  Center.  All  others  should  apply  to  the  National  Technical 
Information  Service. 


If  your  address  has  changed,  or  if  you  wish  to  be  removed  from  the  mailing 
list,  or  if  the  addressee  is  no  longer  employed  by  your  organization,  please 
notify  PL/TSI,  Hanscom  AFB,  MA  01731-5000.  This  will  assist  us  in 
maintaining  a  current  mailing  list. 


Do  not  return  copies  of  this  report  unless  contractual  obligations  or  notices 
on  a  specific  document  requires  that  it  be  returned. 


StCuB'Tv  CA  fiON  OF  This  PAGf 


>•  PGfl  T  S£  CU  «»  T  V  CLA:iSif  iCAT  ION 

Unclasi?if  ied 


:•  seCwRiTv  Classification  AuTnoniTv 


REPORT  DOCUMENTATION  PAGE 


1b  RiSTRiC?  MARK  NOS 


2b  OeCLASStF  iCAT  lON  DOWN G R A  O' NG  SChE  DO UE 


4  PERFORMING  organization  REPORT  NoMSeR^SJ 


3  O'STRiBuT  .ON  A  VAi  lability  QF  REPORT 

Approved  for  public  release; 
distribution  unlimited. 


S  MONlTOP  SG  OnCANiZATlON  .’EPORT  N.VBERiS) 

GL-TR-89-0346 


6«  NAHE  OF  PERFORMING  ORGANIZATION  Bo  OFFICE  SYMBOL  7«  NAME  OF  MONlTOR  NG  ORGANIZATION 

Center  for  Atmospheric  and  i/rapp/.ca6;#> 

Remote  Sounding  Studies  CARSS  Geophysics  Laborator 


6c  ADDRESS  (C|I>.  Slot*  and  Z/P  CiKf*)  Tb  ADDRESS  CiJy.  Slo(»  and  Z/P  Cod»i 

Department  of  Meteorology/CARSS  Hanscom  Air  Force  Base 

University  of  Utah  Bedford,  Massachusetts  01731-5000 

Lake  City,  Utah  84112 


NAME  OF  FUNDING  SPONSOR.NG 
ORGANiZ  AT  lON 


8<  AOwRESS  iCity.  Slate  and  ZIP  C^Jde^ 


80  OFFICE  S'»'M0OC  9  PROCUREMENT  INSTRUMENT  <  DE  F*  T  .  F  .  C  A  T  :  ON  NUMBER 
(If  applik  oblel 

F19628-88-K-0049 


1  title  inciudt  it(ur,iy  ciau  r.cat.nni  Remote  Sounding  of 
ture  Profiles  Usinq  the.. 


ij.  personal  AUTMORISt 
S.  C.  Ou,  K.  N.  Liou 


•3a  type  of  report 

Final 


I3b  Time  covered 


14  DATE  OF  REPORT  i  Yr  .  Mo  .  Da>l  I  15  PAGE  COLNT 

31  December  1989  I  68 


CGSATi  cooes 


18  Subject  terms  fCo^finu#  on  '■nr"*#  if  neceuory  and  identify  by  b<ock  number; 

Remote  Sensing,  Inversion  Theory,  Optical  Measure  Method, 
Differential  Inversion  Method,  HIRS,  Temperature  Retrieval 


'9  A8S  ''RACT  <Cun/inu#  on  i^it'ae  if  necettary  and  identify  b>  b/ocb  number; 

This  report  describes  an  exploratory  investigation  on  the  applicability  of  the 
Optical  Measure  Method  (OMM)  to  the  retrieval  of  atmospheric  temperature  profiles 
using  both  the  simulated  and  HIRS-2  radiances.  We  first  present  the  basic  theory  of 
the  OMM.  Then  we  describe  the  synthetic  retrieval  of  terrperature  profiles  using  HIRS-2 
channels.  Without  performing  an  adjustment  to  the  radiances,  temperature  profiles 
derived  from  both  polynomial  and  polynomial-hyperbolic  functional  fitting  of  radiances 
are  poor.  We  identify  two  major  sources  of  errors  through  a  forward  analysis.  One 
is  the  variation  in  the  sharpness  index  of  the  weighting  function  with  sounding  channels, 
while  the  other  is  the  surface  contribution  to  radiances.  For  these  reasons,  an 
adjustment  scheme  has  been  developed  using  the  concept  of  "scaling  factors",  which 
encompasses  the  effects  of  variation  in  the  sharpness  index,  surface  discontinuity, 
channel  properties,  and  functional  forms.  Retrieved  temperature  profiles  derived  from 


20  OiS'B.BuTlON  A  V  ail  ABILIT  Y  OF  abstract  I 

ai  abstract  security  classification 

jNCL  ASSiF  lED  UNLiMi  TE  O  E  SAME  AS  APT  D  OT  IC  USE  RS  Q 

Unclassified 

.2a  NA./E  OF  RESSO’.S  9i.E  iNC  viDv-Al 

Jean  1.  F.  King 


22t  *€  Lf  --.e  lME  L  R 


22:  G'=  :ES-V9.- 

LY 


UNCLASSIFIED 


;  _  n,Ty  CL  ASS' f  C  ATiON  OF  THIS  »aGE  _  _ 


polynomial-hyperbolic  functional  fitting  of  adjusted  radiances  are  much  inproved.  We  show 
that  the  Differential  Inversion  Method  (DIM)  for  temperature  retrieval  and  OMM  are  math¬ 
ematically  equivalent.  In  practice/  the  DIM  works  better  because  of  the  fitting  of 
radiances  in  the  logarithmic  pressure  scale/  which  accounts  for  some  of  the  surface 
contribution.  Application  of  the  DIM  to  an  archive  of  3473  collocated  temperature  profiles 
and  radiance  data  sets  shows  that  the  accuracy  of  the  retrieved  temperatures  is  within 
about  3  K  below  250  mb.  In  order  to  obtain  a  better  accuracy/  modifications  and  refine¬ 
ments  are  required  for  the  DIM. 

Continuation  of  block  11:  Optical  Measure  Method. 


wis  mki 

fittd  lAB 
tfaannounced 

Justif icatloa 


By - - — j - 

^ptstrlbutlqn/  — 

Availability  Cod<^ 

Bhiail  and/or 


Spaoial 


w 


TABLE  OF  CONTENTS 


Page 

Section  1 

INTRODUCTION  .  1 

Section  2 

OPT I CAL -MEASURE  THEORY  .  3 

2.1  Polynomial  Method  .  3 

2.2  Optical  Number  System  .  5 

2.3  Nonlinear  Hyperbolic  Method  .  7 

Section  3 

APPLICATIONS  OF  THE  OMM  TO  TEMPERATURE  RETRIEVALS  USING 

THE  HIRS  CHANNELS  .  9 

3.1  Determining  of  Weighting  Functions  and 

Radiances  .  9 

3.2  Fitting  of  Radiance  .  11 

3.2.1  Formulation  of  the  PH  Function  .  19 

3.2.2  Prescription  of  mj^  and  m2  .  22 

3.2.3  A  Least-Square  Algorithm  .  22 

3.3  Direct  Application  of  the  OMM  Without  Adjustments  ....  25 

3.4  A  Forward  Analysis  on  the  OMM  .  26 

3.5  Application  of  the  OMM  with  Adjustments  .  34 

3.6  Optimum  Functional  Form  for  the  Planck  Intensity 

Profile  .  38 

Section  4 

COMPARISONS  BETWEEN  THE  DIM  AND  THE  OMM  .  44 

4.1  The  Equivalence  Between  the  DIM  and  the  OMM  .  44 

4.2  Preliminary  Applications  of  the  DIM  and  HIRS  Data  ....  47 

Section  5 

CONCLUSIONS  .  52 

REFERENCES  .  55 

APPENDIX  .  56 

iii 


Section  1 


INTRODUCTION 

There  are  numerous  schemes  that  are  available  for  the  retrieval  of 
temperature  profiles.  In  most  Inversion  approaches,  certain  a  priori  assumptions 
on  the  temperature  profile  or  constraints  on  its  variation  are  required.  The 
root  mean  square  (rms)  errors  of  the  retrieved  temperature  profile  from  the 
traditional  methods  vary  between  1  and  5  K  (Phillips  et  al.  1988).  The  external 
restrictions  generally  take  the  form  of  an  initial  trial  profile  and  the 
departure  allowed  from  the  initial  guess  is  restricted.  The  methods  with  a 
priori  constraint.^  may  work  well  in  an  average  sense  but  large  errors  may  result 
in  extreme  cases. 

The  inference  of  atmospheric  temperature  profile  should  be  based  on  the 
solution  of  the  radiative  transfer  equation.  The  formal  solution  of  the 
radiative  transfer  equation  is  represented  by  the  Fredholm  integral  equation  of 
the  first  kind.  There  are  problems  involving  the  instability  of  tlie  solution  and 
the  ill-conditioned  property  of  the  equation  itself,  since  in  practice  the  data 
available  for  the  performance  of  inversion  are  finite.  Because  the  number  of 
unknowns  are  more  than  the  number  of  equations,  the  solutions  are  not  unique. 

King  (1985)  proposed  a  differential  inversion  method  (DIM),  in  which  the 
Planck  intensity  at  pressure  levels  corresponding  to  the  peaks  of  weighting 
functions  can  be  directly  expressed  in  terms  of  the  sum  of  the  derivatives  of 
upwelling  radiances  in  the  logarithmic  pressure  coordinate.  This  inversion 
technique  that  does  not  involve  a  priori  assumptions  was  developed  based  on  the 
Laplace  transform  of  the  integral  equation.  In  the  follow-up  effort.  King  (1987) 
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developed  a  similar  retrieval  algorithm  referred  to  as  the  optical  measure  method 
(OMM) .  The  fundamental  principle  of  both  the  DIM  and  OMM  is  that  no  initial 
guess  of  the  temperature  profile  or  a  prescription  of  specific  constraints  is 
needed.  The  two  methods  represent  pioneering  efforts  to  break  away  from 
conventional  techniques.  In  the  OMM,  a  certain  functional  form  is  needed  to 
represent  the  Planck  intensity  profile. 

In  our  previous  report,  we  have  explored  the  generalization  and  applicability 
of  the  DIM  by  way  of  a  synthetic  study  (Liou,  et  al.  1989).  We  showed  that  the 
retrieved  temperatures  at  the  lower  four  peak  pressure  levels  of  HIRS  channels 
4-7  are  accurate  to  within  2  K.  In  this  report,  we  present  a  similar 
investigation  on  the  OMM.  In  Section  2,  the  theoretical  foundations  for  the  OMM 
are  presented.  Applications  of  the  OMM  to  synthetic  radiances  of  the  HIRS  15  /xm 
CO2  channels  for  temperature  retrievals  are  described  in  Section  3. 
Subsequently,  the  advantages  and  shortcomings  the  DIM  and  OMM  are  discussed  in 
Section  4.  Finally,  conclusions  are  given  in  Section  5. 
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Section  2 


defined  as 


WjCp/p) 


<ir-(p/p) 


<i<^np) 

and  r-(p/p)  is  the  transmittance  from  the  level  p  to  the  top  of  atmosphere.  We 
substitute  Eqs .  (2.1)  and  (2.2)  into  Eq.  (2.3)  to  obtain 


Z  r,p'  =  Z  b,  p'  r  p'-i  Wj(p)  dp  , 

«=0  f=0  Jo 


(2.5) 


where  p  -  p/p.  The  coefficients  for  the  Planck  function  can  now  be  determined 
by  matching  the  terms  of  equal  power  in  p.  Thus, 


b,  *  r,  £*  p'-i  Wj(p)  dp 


-1 


(2.6) 


and  the  Planck  function  can  be  written  as 


B;(p)  -  Z  r,  [  f  p'-‘  V-(p)  dp 
f=0  I  Jo 


-1 


(2.7'! 


Equation  (2.7)  implies  that  the  Planck  function  can  be  inferred  through 
polynomial  expression  of  p  if  the  coefficients  ,r, ,  are  determinable.  At 
present,  the  number  of  sensing  channels  on  board  satellites  are  limited.  For 
example,  the  HIRS-2  radiometer  has  only  seven  channels  in  the  15  CO^  band  for 
temperature  retrievals.  If  we  use  the  data  from  these  channels,  the  degree  of 
polynomial  expansions  cannot  exceed  six  because  the  coefficients  of  higher- degree 
polynomials  are  not  unique. 

Even  with  all  the  values  of  r,  known,  another  obstacle  inherent  in  the  metliod 
is  the  determination  of  the  integral  in  Eq.  (2.6).  Although  the  integration  from 
p  0  to  p  -  «o  could  be  carried  out  numerically,  its  accuracy  would  be  highly 
questionable.  King  (1987)  proposed  a  new  "optical  number  system"  to  resolve  the 
problem  of  Integration  analytically. 
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2 . 2  Optical  Ntunber  System 

According  to  King  (1985) ,  the  weighting  function  can  be  expressed  in  a 
generalized  form: 


W(P)  =  P  exp(-p  *A) 


(2.8) 


where  k  is  the  "sharpness  index",  and  F  the  Gamma  function.  Two  special  cases 
can  be  derived  from  this  generalized  weighting  function.  For  /c  -  1, 

W(p)  =  p  exp(-p)  (2.9) 

This  expression  corresponds  to  the  weighting  function  derived  from  the  Goody- 
Mayer  random  model.  On  the  other  hand,  for  k  —  2,  we  have 

, -  -  (2.10) 

W(p)  =  ^2/n  p  exp(-pV2) 

This  corresponds  to  the  weighting  function  derived  from  the  Elsasser  band  model. 

Using  the  generalized  weighting  function,  integration  in  Eq.  (2.5)  may  be 
done  analytically.  Using  properties  of  the  Gamma  function,  we  obtain 


r 


p'-^  Wi;(p)dp  =  f(Y7K)  Jo  exp(-pV'')  dP 

x'/*  r((i»l)A] 

° - 


(2.11) 


A  generalized  Gamma  function  may  be  defined  in  the  form 


r  ri(i*i)A] 

’  ■  — ni/^T) — 


(2.12) 


In  this  manner,  Eq .  (2.7)  may  be  expressed  by 


B-(p)  =  S  r,  (r,(f*l)]-i  p' 
^=0 


(2.13) 
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When  <  -  1 ,  the  generalized  Gainina  function  is  reduced  to  the  form  of  Gamma 
function  in  terms  of  natural  number  system.  Since  the  Gamma  function  has  the 


following  property: 


r(i+l)  -  i  r(i) 


for  f  *  1 , 2 . 3 ,  . 


(2.1^) 


we  can  also  construct  an  optical  number  system  n*.  such  that 


r,(iil)  -  i,  r,(i)  ,  for  i  =  1.2,3.  .  . 


(2.15) 


This  optical  number  system  reflects  the  ratio  of  successive  moment  integrals  of 
the  weighting  function,  i.e. 


p'  exp(-p“A)dp 
J"  p*'*  exp(-pV'«)dp 


(2.16) 


For  X  >  1,  which  corresponds  to  a  sharper  weighting  function  with  larger  U- 
and  steeper  slope  on  either  side  of  the  peak,  i  <  i*.  For  x  <  1,  >  2 .  The 

parameter  i,  may  be  referred  to  as  the  sub-natural  and  super-natural  number 
systems  for  i,  <  i  and  i,  >  i,  respectively.  For  retrieval  purposes,  the  sub- 
natural  number  system  Is  desired,  because  sharper  weighting  functions  would  giv^ 
more  accurate  retrieval  results  (Liou  and  Ou,  1989). 

In  practice,  there  is  a  weighting  function  profile  for  each  sounding  channel. 
A  distinctive  value  of  x  can  be  assigned  to  each  channel  by  fitting  the 
generalized  form  to  the  weighting  function.  Thus,  if  we  evaluate  B-(p)  according 
to  Eq .  (2.7),  we  will  obtain  several  values  of  B-(p)  that  correspond  to  different 
X.  A  conversion  of  B-(p)  to  temperature  may  not  lead  to  a  uniform  x  value  at  the 
same  pressure  level.  To  determine  a  mean  temperature  profile,  a  weighting  method 
is  usvially  used  (Smith,  19/0). 
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2 . 3  Nonlinear  Hyperbolic  Method 

The  key  to  the  success  of  the  OMM  lies  in  the  adequate  fitting  of  a  finite 
number  of  radiance  measurements  by  a  smooth  function  of  a  combination  of  smooth 
functions.  Higher  -  degree  polynomials  may  fit  every  radiance  data  point  v.-ell,  but 
thc'y  may  also  carry  spurious  oscillations  between  data  points.  This  will  lead 
to  unrealistic  solutions  of  the  Planck  function.  Thus,  we  are  forced  to  seek  tk.e 
fitting  of  radiances  that  will  not  produce  unnecessary  wiggles  between  data 
points,  but  at  the  same  time  the  fitted  curve  will  match  each  data  point  •..i'li 
reasonable  accuracy. 

One  possible  candidate  for  the  radiance  fitting  is  the  nonlinear  liyporbolic 
function.  According  to  Chandrasekhar  (1950),  the  upwelling  intensity  iro:r.  a 
plane-parallel  atmosphere  can  be  expressed  as  a  sum  of  hyperbolic  functions 
Tlius,  R[p(i/)]  may  be  expressed  in  the  form 


R[p(i.)l  =  S  — ^ 
j=0  l^nijp 


(2,17) 


.Actually,  Eq .  (2.17)  is  equivalent  to  a  combination  of  geometrical  series: 


.2  =  .S  -m,p  *  (m,p)^  - 

j=0  l+m^p  j=0 


(2 . 18^ 


Comparing  Eq .  (2.18)  with  Eq .  (2.1),  we  note  that 


r,  =  S  r,  (-nij)' 
j=0  '' 


(2.19) 


and  the  coefficients  of  the  Planck  function  are  given  by 


b,  = 


’=0 


(2 . ?0) 
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It  follows  that  Eq.  (2.21)  may  be  generalized  in  the  form 


B;(p)  *  2  r.  exp^(-m,p) 

j=o  •’ 


(2.24) 


Using  the  combination  of  hyperbolic  expressions  for  radiances,  we  obtain  a 
series  expansion  of  B-  in  terms  of  the  optical  exponential  function.  Evaluation 
of  the  optical  exponential  function  must  be  done  through  Eq .  (2.23).  For  x  -  1 , 
this  equation  reduces  to  the  exponential  function. 
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Section  3 


APPLICATIONS  OF  THE  OMM  TO  TEMPERATURE  RETRIEVALS 
USING  THE  HIRS  CHANNELS 


In  this  section,  we  shall  present  the  results  of  applying  the  OMM  to 
temperature  retrievals  by  performing  synthetic  analyses  using  the  properties  of 
HIRS  channels  in  the  15  /im  CO2  band. 


1 . 1  Determination  of  Weighting  Functions  and  Radiances 

There  are  seven  channels  in  the  HIRS  15  (im  CO2  band  with  central  wavenumbers 
at  668,  679,  690,  702,  716,  732,  and  748  cm'^.  The  smaller  wavenumber 
corresponds  to  broader  weighting  functions  with  maxima  located  at  higher 
altitudes.  The  last  two  channels  are  also  affected  by  the  absorption/emission 
due  to  water  vapor.  Details  of  the  channel  characteristics  are  listed  in  Table 
1  . 

Table  1.  HIRS  Channel  Characteristics 


Channel 

1/  (cm'^) 

Ai/ 

Principal 

Absorbers 

Level  of 
U 

1 

668 

666 

670 

4 

CO2 

30 

2 

679 

674 

684 

10 

CO2 

60 

3 

690 

685 

697 

12 

CO2 

100 

4 

702 

696 

712 

16 

CO2 

250 

5 

716 

708 

724 

16 

CO2 

500 

6 

732 

724 

740 

16 

CO2/H2O 

750 

7 

748 

740 

756 

16 

CO2/H2O 

900 
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We  used  the  CO2  absorption  coefficients  computed  by  Chou  and  Kouvaris  (1986) 
based  on  line-by-line  data  compiled  by  Rothman  et  al,  (1983)  to  compute 
transmittances  and  weighting  functions,  using  the  k- correlated  method.  This 
method  transforms  the  spectral  integral  over  the  wavenumber  space  into  an 
integral  over  the  domain  of  the  accumulated  frequency  distribution  of  the 
absorption  coefficient.  Furthermore,  to  account  for  the  inhomogeneous 
atmosphere,  we  assume  that  the  variation  in  the  absorption  coefficient  depends 
only  on  pressure  and  temperature  so  that  the  frequency  distributions  of  the 
absorption  coefficient  can  be  correlated  through  pressure  and  temperature.  In 
practice,  we  use  the  cumulated  frequency  distribution  (g)  look-up  table  to 
determine  the  relationship  between  the  absorption  coefficient  at  an  arbitrary 
level  and  that  at  the  reference  level.  The  reference  absorption  coefficient  is 
determined  from  each  g  value.  The  spectral  transmittance  is  obtained  by 
numerical  integration  over  the  g  space  (Liou  and  Ou,  1988). 

A  numerical  method  was  developed  to  fit  the  HIRS  weighting  functions  to  the 
generalized  form  denoted  in  Eq.  (2.8)  to  obtain  the  sharpness  index,  k, 
associated  with  each  channel.  It  is  similar  to  the  least-square  method. 
However,  the  quantity  that  is  to  be  minimized  is  the  weighted  square  sum  of 
errors.  More  weight  is  placed  on  the  error  near  the  peak  of  the  weighting 
function,  because  a  significant  portion  of  the  upwelling  radiance  comes  from 
emission  near  the  peak  level.  The  resulting  nonlinear  algebraic  equation  is 
solved  by  Newton's  iteration  method.  The  sharpness  index  of  the  fitted  weighting 
functions  for  Channels  1-7  are  0.49,  1.56,  1.50,  2.19,  2.34,  4.34,  and  3.16 
respectively.  Except  channel  1,  the  k  values  associated  with  all  other  channels 
are  larger  than  1.  The  k  value  deviates  from  1  and  2,  implying  that  neither  the 
Goody  random  band  model  nor  the  Elsasser  regular  band  model  can  properly  simulate 
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the  total  transmittance  in  the  15  /im  CO2  band.  The  errors  for  the  fitted 
weighting  functions  are  within  +  0.02.  Since  the  fitting  method  places  more 
weight  on  minimizing  errors  near  the  peaks  than  in  the  wings,  the  percentage 
errors  near  the  peaks  (<  2%)  is  much  smaller  than  those  in  the  wing  regions. 
Figures  1(a) -(g)  show  the  weighting  functions  computed  from  the  k-correiattd 
method  and  the  fitted  generalized  weighting  functions.  The  fitting  for  Channel 

1  is  not  as  satisfactory  as  those  for  other  channels,  due  to  the  secondary  peak 
above  10  mb.  The  figures  also  show  that  the  weighting  function  of  Channel  1  is 
much  broader  than  other  channels.  Channels  6  and  7  are  the  sharpest  channels  with 
peaks  close  to  the  surface.  The  fitted  weighting  functions  for  these  two  channels 
are  slightly  broader  than  those  computed  from  the  k-correlated  method.  The 
present  results  were  obtained  by  the  systematic  optimization  method,  and  therefore 
can  be  applied  to  any  weighting  function. 

The  channel  radiance  values  were  then  obtained  from  synthetic  computations. 
We  used  the  US  Standard  Atmospheric  Temperature  Profile  and  evaluate  the  radiances 
according  to  Eq.  (2.3).  The  generalized  weighting  functions  were  used.  Table  2 
lists  the  computed  upwelling  radiances  for  the  first  seven  HIRS  channels.  Wo 
notice  that  these  radiances  increase  from  the  band  center  (Channels  2  and  3)  to 
the  band  wing  (Channels  6  and  7). 

3 . 2  Fitting  of  Radiance 

The  seven  upwelling  radiances  were  fitted  to  either  a  5th-degree  polynomial 
function  or  a  combined  polynomial -hyperbolic  (PH)  function.  In  the  case  of  the 
5th-degree  polynomial  fitting,  we  used  the  subroutine  DRCURV  in  the  Internationa] 
Mathematical  and  Statistical  Library  (IMSL).  The  resulting  curve  is  shown  in  Fig. 

2  as  the  solid  line.  There  are  two  reasons  for  choosing  the  5th-degree 
polynomial.  For  a  lower-degree  polynomial,  the  fitting  will  not  be  exact.  On  the 
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Fig.  1(a).  The  weighting  function  of  HIRS-2  Channel  1.  The  solid  and 

dashed  lines  represent  results  computed  from  the  correlated  k- 
distribution  method  and  the  generalized  weighting  function, 
respectively . 


1(d). 


Same  as  Fig. 


except  for  HIRS-2  Channel 


.2  0.3  0 

NEXOHTING  FUNCTION 


Fig.  1(e).  Same  as  Fig.  1(a),  except  for  HIRS-2  Channel  5 
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PRESSURE. 


Fig.  Kg).  Same  as  Fig.  1(a),  except  for  HIRS-2  Channel  7. 


Table  2.  Radiance  values  for  HIRS  channels. 


Channel 

Radiance 

(erg/cm^-cm'^-  sec  -sr) 

1 

53.5 

2 

46.0 

3 

44.4 

4 

54.9 

5 

60.1 

6 

70.6 

7 

78.2 

other  hand,  for  a  higher-degree  polynomial,  there  will  be  iinrealiEtic 
oscillations  between  data  points.  Figure  2  shows  that  the  fitted  curve  has  a 
local  minimum  between  p  -  100  and  200  mb.  The  curve  fits  almost  exactly  to  all 
the  radiances  except  those  of  Channels  2  and  3.  In  addition,  there  are  slight 
but  noticeable  oscillations  between  Channels  4  and  7. 

In  the  case  of  fitting  the  radiances  to  the  PH  function,  we  have  developed 
;in  algorithm  based  on  the  principle  of  the  least-square  method.  This  is 
described  in  the  following  subsections. 

3.2.1  Formulation  of  the  PH  Function 

The  combined  PH  function  can  be  written  as  follows: 

-  _  _  Ni  _  N2  L 

R(p(«^)l=  S  r,p'-‘  *  Z  —!-=  ,  ,  ,  p 

i  =  l  i  =  l  l+ro,p 

which  is  a  combination  of  Eq.  (2.1)  and  (2.17).  The  selection  of  such  a  form  was 
based  on  reasons  related  to  the  distribution  of  radiances.  Since  the  peak  of  the 
Channel  1  weighting  function  is  well  within  the  stratosphere,  in  which  the 
temperature  is  generally  near  isothermal  or  increases  with  altitude,  the  measured 
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polyriMiial 
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fitting  of  a  5th-order  polynomial  (solid  line)  and  a  polynomial- 
hyperbolic  function  (dashed  line)  to  the  synthetic  radiances  for 
the  seven  HIRS-2  channels.  The  US  Standard  Atmosphere  Temperature 
profile  was  used  in  s)mthetic  computations. 


Channel  1  radiance  is  usually  higher  than  Channel  2  radiance.  The  radiances 
increase  from  Channel  3  to  Channel  7  as  the  peaks  of  the  weighting  function 
gradually  lower  toward  surface.  The  corresponding  temperature  at  the  peak 
increases.  Usually,  radiances  measured  by  Channel  2  or  Channel  3  are  the 
smallest.  This  implies  that  the  fitted  curve  must  be  able  to  generate  such  local 
minima.  Thus,  a  second-  or  higher-degree  polynomial  is  required.  However,  a' 
the  degree  of  the  fitted  polynomial  increases,  so  does  the  amplitude-  of 
fluctuation  between  data  point,  which  may  cause  serious  instability  on  the 
solution.  An  alternative  is  to  fit  the  radiances  to  a  series  of  non- linear- 
hyperbolic  functions  denoted  in  Eq .  (2.17).  However,  our  preliminary  studies 
show  that  the  resultant  fittings  based  on  the  least-square  principle  are  far  from 
satisfactory.  One  serious  defect  is  that  it  cannot  reproduce  the  local  minimum 
of  Channel  2  radiances,  because  the  functional  form  1/(1  +  x)  is  monotonous  on 
either  side  of  x  -  -1. 

Finally,  we  find  that  the  combined  PH  function  can  simulate  the  local 
minimum,  and  at  the  same  time  it  meets  the  smoothness  requirement.  Since  the 
number  of  data  points  in  the  present  case  is  seven,  the  fitting  should  have  no 
more  than  seven  undetermined  coefficients.  If  the  number  of  unknown  coefficien.ts 
exceeds  seven,  the  fitting  is  not  unique  unless  additional  constraints  are 
imposed.  Thus,  for  best  fitting,  we  can  have  only  three  combinations:  (1)  N, 
-  1,  N2  -  3  (7  unknowns),  (2)  Nj  -  Nj  -  2  (6  unknowns),  and  (3)  Ni  -  3,  N2  -  2  (7 
unknowns).  As  stated  above,  case  (1)  is  the  sum  of  pure  hyperbolic  functions, 
which  cannot  reproduce  the  local  minimtim.  At  present,  we  have  selected  case  (2) 
IS  a  preliminary  test  of  the  functional  form  since  the  coefficients  in  this  c-jst 
I'an  be  easily  determined. 
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3.2.2  Prescription  of  aH]  and 


From  Eq.  (3.1),  let  Nj  -  Nj  -  2.  We  have 


R[p(i/)]  =  ri  ♦  rjP  + 


1  +  m^p  1  +  mjP 


(3.2) 


where  r^,  V2,  L^,  L2,  mi,  and  m2  are  coefficients  to  be  determined. 

Equation  (3.2)  Is  a  non-linear  function,  if  m^  and  m2  are  considered  to  be 
unknown  coefficients.  However,  if  we  prescribe  mj  and  m2,  then  it  becomes  a 
linear  combination  of  three  functional  forms  with  one  constant.  The  remaining;, 
coefficients  can  be  found  from  a  least-square  method.  In  general,  m^  and  m2  can 
be  real  numbers  from  -<»  to  +«.  Fortunately,  by  carefully  reviewing  the 
derivation  of  the  retrieval  formula,  we  find  that  in  order  for  the  method  to  be 
mathematically  consistent,  mj  and  m2  must  be  within  a  certain  range.  We  notice 
that  the  equality. 


—  *l-x*x2-x^...  =  E  (-X)' 

1  +  X  £=0 


(3.3) 


is  valid  only  when  |x|  <  1.  Applying  this  restriction  to  Eq.  (2.18),  we  conclucio 
that  |mip|  and  lm2p|  must  be  less  than  1,  or  |mi|  and  |m2l  must  be  less  than  1/p. 
Otherwise,  Eqs.  (2.18)  and  (2.24)  will  not  be  valid.  Since  p  ranges  between  0 
and  1000  mb,  |mi|  and  |m2|  have  to  be  less  than  10'^.  This  small  range  allows  us 
to  pre-select  a  few  values  of  mj  and  m2  within  this  range  to  compute  the  best- fit 
coefficients  in  Eq.  (3.2)  and  the  corresponding  rms  error. 


3.2.3  A  Least-Square  Algorithm 

Based  on  Eq.  (3.2),  the  square  sura  of  the  difference  between  the  fitted  and 
the  actual  values  of  R(p(i/)]  may  be  expressed  by 
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(3.4) 


E  =  S  (rj  ^  rj  Pj  +  - -  Rj)2 

j  =1  1  4-  nil  pj  1  4.  rn2  Pj 

To  minimize  E,  we  set  the  partial  differentials  with  respect  to  the  coefficients 
to  zero ,  vis . , 


dE  _  dE  ^  dE  _  dE 
'5x1  ~  ^  ^  ^ 


0 


(3.5; 


This  minimization  process  leads  to  a  system  of  four  linear  equations,  which  can 
be  expressed  in  the  following  matrix  form: 


where 


7  P  Xi  Xj 

T  p2  x^ 

x^p  xl  x^ 

T2X2P  ^2  xi 


Xi  =  .2  1/(1  +  niipp 

J=1 


ri 

R 

^2 

Li 

Riq 

L2 

W2 

(3.6) 


_  7 


Xf  =  .2  1/(1  4 
J=1 


1  =  1.2 


1  =  1.2  , 


t  _ 

X^  =  1/[(1  4  miPj)  (1  4  ra2Pj)] 


7  _ 
P  =  S  Pj 
j»l 


p2  =  S  Pj 

j=l 
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_  7  _ 

XjP  =  pj/(l  +  m^pj).  i  =  1,  2  , 

7 

R  =  E  Rj  , 

j=l  ' 

_  7  _ 

7 

RXi  =  ^^1  *  “iPj>  •  i  =  1.  2  . 

Solution  of  ri,  r2,  L^,  and  L2  depends  on  the  prescribed  value  of  mj  and  m2.  The 
rms  error  is  defined  by 

e  =  ^/e/7 

Since  the  values  of  mj  and  m2  should  be  between  -10'®  and  10'^  ,  we  select  six 
representative  values:  ±  3.5  x  lO"'",  ±  6.5  x  10'*,  and  ±  9.5  x  10'*.  For  each 
pair  of  mi  and  m2,  we  may  solve  the  linear  system  of  Eq.  (3.6)  and  compute  the 
rms  error  from  Eq.  (3.7).  The  idea  is  to  determine  which  pair  of  (mi, m2)  produce 
minimum  rms  error.  Table  3  lists  the  rms  error  for  all  pairs  of  (mi, m2)  using 
the  synthetic  radiances  listed  in  Table  2.  Clearly,  Table  3  shows  that  the  rms 
error  is  generally  smaller  for  mira2  <  0  (one  negative)  than  for  mim2  >  0  (both 
positive  or  negative).  Furthermore,  for  the  case  mim2  <  0,  the  rms  error  for  the 
pair  of  (mi, m2)  is  the  same  as  for  the  pair  of  (m2, mi).  The  minimum  rms  error  is 
2.8  for  mi  -  9.5  x  10  *  and  m2  -  -3.5  x  10'*,  and  the  coefficients  in  this  case 
are  ri  -  79.6,  r2  -  0.A98,  Li  -  AA5  and  L2  —  -473.  The  fitted  PH  function  is  also 
shown  in  Fig.  2  as  the  dashed  line.  Although  the  general  fitting  of  the  PH 
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Table  3.  The  root -mean- square  error  of  radiances  for  selected  pairs  of  (mi, m2) 


m2  X  10* 

-9.5 

-6.5 

mj  X  10* 

-3.5  3.5 

6.5 

9.5 

-9.5 

4.2 

3.8 

3.7 

3.5 

3.4 

3.4 

-6.5 

3.8 

4.1 

3.5 

3.3 

3.2 

3.1 

-3.5 

3.7 

3.6 

3.9 

3.2 

2.9 

2.8 

3.5 

3.5 

3.3 

3.2 

3.7 

23.4 

13.9 

6.5 

3.4 

3.2 

2.9 

8.5 

3.7 

4.6 

9.5 

3.4 

3.1 

2.8 

5.1 

3.2 

3.6 

function  is  less  satisfactory  than  the  polynomial  function,  particularly  near  the 
local  minimum,  the  PH  function  is  smoother  than  the  polynomial  function.  Better 
fittings  can  be  obtained  with  |m,  |  >1  (Leon  and  King,  1989).  However,  as  noted 
above,  the  series  will  not  converge  and  the  OMM  is  not  mathematically  valid  in 
this  case. 

3 . 3  Direct  Application  of  OMM  Without  Adjustments 

Retrieval  exercises  were  carried  out  based  on  the  two  functional  forms  for 
radiances.  Planck  functions  derived  from  the  polynomial  form  can  be  written  as 


8(1/1, p) 


6  r,  p'-i 
ill  T7^ 


(3.8) 


where  k  is  associated  with  the  ith  channel,  Planck  functions  based  on  the  PH 
form  can  be  written  as  a  combination  of  polynomial  and  optical -exponential 
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functions  as  follows: 


B(i/i,p)  =  Ti  +  rjp/l*!  +  Li  exp,(-mip)  +  Lj  exp^(-ni2P) 

In  this  way,  since  there  are  seven  channels,  each  with  a  different  value  for  k, 
we  may  obtain  seven  profiles  for  B(I/£,p).  Each  Planck  profile  can  be  converted 
to  a  temperature  profile  T(Vi,p).  We  used  the  method  of  weighted  average  (Smith, 
1970)  to  obtain  the  final  temperature  profile  as  follows: 

-  7  _  7  _ 

T(p)  =  S  [W(v,,p)  T(Vi,p)]  /  2  W(v,,p)  . 

i=l  i=*l 


Figure  3  shows  the  retrieved  temperature  profiles  based  on  Eqs.  (3.8)  and 
(3.9),  along  with  the  reference  standard  temperature  profile.  Both  retrieved 
profiles  deviate  significantly  from  the  standard  temperature  profile.  The 
profile  derived  from  the  PH  function  is  better  than  that  from  the  polynomial. 
A  further  examination  of  the  figure  reveals  that  the  profile  from  the  PH  function 
is  colder  by  5-12  K  below  250  mb  but  warmer  by  0-5  K  above  this  pressure  level. 

In  an  attempt  to  investigate  the  reason  for  discrepancies,  we  found  two 
problems  that  need  to  be  resolved  in  order  for  the  OMM  to  be  practical  for 
temperature  retrievals.  The  first  problem  is  related  to  the  variation  in  k. 
while  the  second  is  associated  with  the  discontinuity  of  the  temperature  profile. 
These  problems  were  identified  through  a  forward-analysis  procedure  which  will 
be  described  in  the  following  subsection. 


3 . 4  A  Forward  Analysis  on  the  OMM 

In  this  section,  a  forward  analysis  will  be  performed  to  illustrate  the 
reasons  for  the  deviation  of  the  retrieved  temperature  profile  from  the  standard 
temperature  profile. 
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3.  Retrieved  temperature  profiles  based  on  the  fitting  of  synthe 
radiances  shown  In  Fig.  2. 


According  to  the  description  in  Section  2,  the  Planck  profile  can  be 
represented  by  either  polynomial  functions  or  optical  exponential  functions  or 
a  combination  of  the  two  as  shown  by  Eq.  (3.9).  It  is  not  known  how  well  such 
functional  forms  can  match  the  Planck  profiles.  In  the  previous  subsection,  we 
demonstrated  that  the  PH  function  yields  better  retrieval  results  than  fix' 
polynomial  function.  For  this  reason  the  PE  function  defined  by  Eq.  (3.9)  i:; 
chosen  to  fit  the  standard  Plank  profile. 

Given  that  m^  and  ra2  are  prescribed,  we  use  the  least-square  method  similar 
to  the  one  described  in  Sec.  3.2.3  to  perform  numerical  fittings.  We  first 
define  a  square  sum  of  errors  as 


N 

E  *  S  [r^  +  rj  Pj/1,!  +  Lj  exp„(-miPj) 

*  Lj  exp„(-m2Pj)  -  Bj]2  . 


(3.11) 


Performing  the  minimization  procedure,  we  obtain  the  following  linear  system: 


■ 

■ 

. 

N  n  xi  xz 

n  xiH  Xz^ 

^1 

rz 

B 

Bn 

Xi  ATin  Xi  XiXz 

Li 

Wi 

Xz  XzH  XiXz  xi 

Lz 

where 


_  ^ 

Xi  »  2  exp,,(-miPj)  ,  i  =  1,  2 

j=l 


N 

^  *  S  exp„(-miPj)5 


1  »  1,  2 


XiXz 


N 

S 


j=l 


exp„(-mjPj)  exp„(-m2Pj) 
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N 

n  =  pj/ij 


N 

2  (Pj/iJ)" 

j=i 


_  N 

Xi^  “  Pj  exp,(-mjpj) 


i  *  1,  2 


N 

B  =  S  B.  , 

j=l 

N 

Bn  =  .S  Bjpj/l^!  , 

J=1 

_  N 

*  Z  Bj  exp^(-iniPj)  ,  i  »  1,  2  , 

j  »1 

In  computing  the  preceding  matrix  elements,  the  evaluation  of  the  optical 
exponential  function  was  done  by  the  series  summation  method  according  to  Eq. 
(2.23).  We  found  that  about  15  terms  are  sufficient  for  the  series  to  converge 
to  a  constant  number.  Since  each  channel  has  a  different  Planck  intensity 
profile,  and  since  the  optical  exponential  function  varies  with  the  sharpness 
index  k,  the  coefficients  r^,  r2,  L^,  and  L2  depend  on  both  the  spectral  band  of 
each  channel  for  k.  For  the  seven  HIRS  channels,  there  are  49  sets  of 
coefficients.  These  coefficients  are  listed  in  Table  4.  It  is  evident  that  for 
the  same  k,  the  coefficients  vary  slightly  with  channels.  However,  for  the  same 
channel,  the  coefficients  depend  strongly  on  k.  Also,  these  sets  of  values  are 
different  from  the  values  in  Section  3.2.3.  The  first  problem  can  now  be  stated 
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Table  4.  Coefficients  of  the  fitted  polynomial -hyperbolic  functional  form  for  the 
seven  HIRS  channels  in  terms  of  different  sharpness  indices,  k 


Channel 

1C 

0.49 

1.56 

1.50 

2.19 

2.34 

4.34 

3.16 

0.315E+04 

0.131E+04 

0.186E+04 

0.632E+03 

0.575E+03 

0.187E+03 

0.336E+03 

1 

0.250E+01 

0.146E+01 

0.195E+01 

0.933E+00 

0.894E+00 

0.537E+00 

0.690E+00 

"^1 

0.989E+03 

0.818E+03 

0.107E+04 

0.598E+03 

0.585E+03 

0.413E+03 

0.493E+03 

b 

-0.409E+04 

-0.208E+04 

-0.288E+04 

-0.118E+04 

-0.112E+04 

-0.554E+03 

-0.783E+03 

0.306E+04 

0.130E+04 

0 . 184E+04 

0.626E+03 

0.560E+03 

0.183E+03 

0.331E+03 

2 

0.236E+01 

0.145E+01 

0.154E+01 

0.928E+00 

0.890E+00 

0.534E+00 

O.886E1OO 

0.911E+03 

0.816E+03 

0.107E+04 

0.598E+03 

0.583E+03 

0.411E+03 

0.491Et03 

^2 

-0.393E+04 

-0.207E+04 

-0.286E+04 

-0.118E+04 

-O.lliE+04 

-0. 548E+03 

-0.778E403 

0.302E+04 

0.128E+04 

0.182E+04 

0.617E+03 

0.560E+03 

0.179E+03 

0.325EI03 

3 

0.236E+01 

0. 144E+01 

0.193E+01 

0.923E+00 

0.883E+00 

0.531E+00 

0.681EH00 

^1 

0.950E+03 

0.813E+03 

0.106E+04 

0.593E+03 

0.580E+03 

0.410E+03 

0.489Ef03 

■-2 

-0.409E+04 

-0.205E+04 

-0.285E+04 

-0.117E+04 

-O.llOE+04 

-0.545E+03 

-0.770E+03 

0.287E+04 

0.127E+04 

0.181E+04 

0.606E+03 

0.549E+03 

0.173E+03 

0.319E+03 

U 

0.240E+01 

0.143E+01 

0.191E+01 

0.913E+00 

0.874E+00 

0.526E+00 

0.676E+00 

"-I 

0.973E+03 

0.808E+03 

0.105E+04 

0.588E+03 

0.576E+03 

0.407E+03 

0.486E+03 

■-2 

-0.405E+04 

-0.203E+04 

-0.281E+04 

-0.115E+04 

-0.108E+04 

-0.539E+03 

-0.762E403 

0.287E+04 

0.125E+04 

0.178E+04 

0.596E+03 

0.540E+03 

0. 168E+03 

0.312E+03 

5 

0.218E+01 

0.141E+01 

0. 189E+01 

0.904E+00 

0.866E400 

0.521E+00 

0.669E400 

"'I 

0.957E+03 

0.799E+03 

0.105E+04 

0.584E+03 

0.572E+03 

0.404E403 

0.482F.  +  03 

4 

-0.368E+04 

-0.201E+04 

-0.279E+04 

-0.114E+04 

-0.107E+04 

-0.532E+03 

-0.753Et03 

0.287E+04 

0. 122E+04 

0.175E+04 

0.581E+03 

0.528E+03 

0.162E+03 

0.302E+03 

6 

0.227E+01 

0.139E+01 

0.186E+01 

0.891E+00 

0.854E+00 

0.513E+00 

0.658E400 

L? 

0.970E+03 

0.790+E03 

0.103E+04 

0.577E+03 

0.566E+03 

0.399E+03 

0.477E+03 

‘•2 

-0.362E+04 

-0.197E+04 

-0.274E+04 

-0.112E+04 

-0.105E+04 

-0.522E+03 

-0.740E+03 

I"! 

0.271E+04 

0. 120E+04 

0. 171E+04 

0.568E+03 

0.514E+03 

0. 135E+03 

0.293E+03 

7 

0.223E+01 

0. 138E+01 

0.184E+01 

0.879E+00 

0.841E+00 

0.506E+00 

0.650E+00 

0.984E403 

0.784E+03 

0.102E+04 

0.571E+03 

0.558E+03 

0.395E+03 

0.471Ef03 

*•2 

-0.363E+04 

-0.195E+04 

-0.270E+04 

-O.llOE+04 

-0.103E+04 

-0.512E+03 

-0.727E103 
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as  follows:  "Due  to  the  variation  in  the  sharpness  index  with  the  channel,  it 
appears  not  practical  to  use  only  one  set  of  coefficients  to  represent  all  the 
measured  radiances . " 

Since  each  radiance  value  corresponds  to  a  different  k,  using  a  single  set 
of  coefficients  to  represent  all  the  radiances  as  has  been  shown  in  Subsection 
3.2.3  would  produce  large  temperature  errors.  In  the  retrieval  exercise 
described  in  Subsection  3.3,  none  of  the  retrieved  temperature  profiles  is  close 
to  the  standard  profile. 

To  show  how  well  the  PH  function  defined  by  Eq.  (3.11)  fits  the  input  Planck 
intensity  profile,  the  rms  errors  for  the  Planck  function  and  temperature  are 
presented  in  Table  5.  Also  in  this  table  are  temperature  rms  errors.  For  x  < 
1,  the  fitting  is  poor  because  the  matrix  approaches  singular.  On  the  other 
hand,  for  k  >  2,  the  temperature  rms  errors  are  on  the  order  of  1.2  -  2.0  K. 
This  indicates  that  the  PE  function  is  potentially  useful  to  produce  real 
atmospheric  Planck  profiles. 

Next,  we  examine  the  errors  in  the  radiances  generated  by  the  fitted  Planck 
profile.  We  compute  three  sets  of  "radiances".  The  first  is  based  on  the  input 
Planck  profile  with  a  contribution  from  the  surface  separate  from  the  integral 
as  follows; 

Ri(Pi)  =  BpCp,)  rp(p/Pi) 

+  IJ*  Bp(p)  W-(k^,p/p^)  dp/p  (3.13) 

Note  that  Eq .  (3.13)  is  equivalent  to  Eq.  (2.3),  if  the  temperature  for  p  >  Ps 
is  assumed  to  be  isothermal.  The  second  set  of  radiances  was  obtained  from  the 
following  equation: 

Ri(Pi)  =  Bp(p,)  rp(p/Pi) 

+  £**  Bp(p)  W-(fCj,p/Pi)  dp/p  (3.1^) 
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Table  5. 

The  rms  error  for 
channels  in  terms 

the  fitting  of  Planck  profiles  of  the 
of  different  sharpness  indices,  k 

first 

seven  HIRS 

Channel 

K 

0.49 

1.56 

1.50 

2.19 

CM 

1 

4.34 

3.16 

1  Sb 

19.1 

1.9 

4.1 

1.5 

1.5 

1.5 

1.5 

14.6 

1.7 

3.5 

1.4 

1.3 

1.3 

1.3 

2  Tb 

9.5 

1.7 

2.1 

2.1 

1.5 

1.5 

1.5 

AT 

7.4 

1.5 

1.7 

2.0 

1.4 

1.3 

1.3 

3  AB 

** 

2.7 

6.2 

1.4 

1.6 

1.4 

1.5 

** 

2.4 

5.6 

1.3 

1.5 

1.3 

1.3 

4  AB 

** 

1.7 

2.9 

1.5 

1.6 

1.4 

1.4 

AT 

** 

1.7 

2.4 

1.4 

1.5 

1.3 

1.3 

5  AB 

** 

1.5 

1.7 

1.4 

1.4 

1.4 

1.4 

** 

1.4 

1.3 

1.4 

1.4 

1.3 

1.3 

6  AB 

** 

2.0 

2.4 

1.4 

1.5 

1.4 

1.4 

1.9 

2.2 

1.3 

1.5 

1.3 

1.2 

7  AB 

21.0 

1.9 

4.1 

1.6 

1.4 

1.3 

1.3 

AT 

20.5 

1.9 

4.0 

1.5 

1.3 

1.2 

1.2 

**  Error  too  big  due 

to 

unstable  matrix 

solution 

32 


where  B-(p)  Is  the  fitted  Planck  intensity  profile.  The  third  set  of  radiances, 
^3(Pi) .  was  obtained  from  Eq.  (3.2)  using  the  coefficients  obtained  in  the 
forward  analysis .  The  "surface  contribution"  in  this  case  is  obtained  from  the 
difference  between  R3(Pi)  and  the  atmospheric  contribution  to  the  radiance  as 
follows : 

RsCPi)  (Surface  Contribution) 

=  R3<Pi)  -  B-(p)  W-(«j,p/pj)  dp/p  (3.15) 

A  comparison  of  Ri(pi)  and  R2(Pi)  will  reveal  the  effect  of  errors  in  B-(p)  for 
p  <  Pj,  while  the  difference  between  R2(Pi)  and  R3(Pi)  will  show  the  impact  of  the 
irregular  behavior  of  B-(p)  for  p  >  Ps-  Table  6  lists  the  radiances  that  have 
been  computed  for  the  first  seven  HIRS  channels,  where  the  numbers  in  parentheses 
are  associated  with  surface  contributions.  The  errors  in  radiances  produced  by 
B-(p)  for  p  <  p»  are  generally  less  than  5  erg/cmVsec/cm* Vsr .  The  larger 
differences  for  the  first  three  channels  are  due  to  the  fitting  errors  in  B-(p) 
in  the  stratosphere.  Surface  contributions  are  significant  for  the  last  three 
channels  leading  to  different  R2(Pi)  and  R3(pi)  .  This  is  because  B-(p)  behaves 
irregularly  for  p  >  pg,  where  the  values  of  weighting  functions  are  still 
comparable  to  the  peak  value.  Clearly,  errors  will  be  generated  in  the 
temperature  retrieval  due  to  the  limitation  of  the  functional  form  associated 
with  B-(p) .  Thus  the  second  problem  can  be  stated  as  follows:  "Due  to  the 
irregular  behavior  of  the  fitted  functional  form  of  B-(p) ,  retrieved  temperatures 
near  the  surface  and  in  the  stratosphere  will  suffer  large  errors. 

To  make  the  OMM  a  feasible  scheme  for  temperature  retrievals,  certain 
modifications  or  adjustments  are  required  in  order  to  circumvent  the  preceding 
problems . 
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Table  6.  The  radiances  for  the  first  seven  HIRS  channels,  where  the  numbers 


in  parentheses 
of  radiances  is 

are  the 
erg/cm^ 

respective 

/sec/cm‘^/sr 

surface 

contributions.  The 

uti  i  1 

Channel 

1 

2 

3 

4 

5 

6 

7 

K 

0.49 

1.56 

1.50 

2.19 

2.34 

4.34 

3.16 

53.5 

46.0 

44.4 

54.9 

60.1 

70.6 

78.2 

(0.0) 

(0.0) 

(O.C) 

(0.0) 

(0.2) 

(9.6) 

(28.4) 

47.6 

42.7 

37.6 

53.8 

60.0 

70.3 

77.4 

(0.0) 

(0.0) 

(0.0) 

(0.0) 

(0.2) 

(9.4) 

(27.5) 

48.2 

42.7 

37.6 

54.3 

64.7 

75.1 

67.0 

(0.0) 

(0.0) 

(0.0) 

(0.0) 

(4.9) 

(14.2) 

(17.1) 

3 . 5  Application  of  OMM  with  Adjustments 

Based  on  the  forward  analysis  described  in  the  previous  subsection,  we 
conclude  that  the  problems  of  the  variation  in  the  sharpness  index  of  the  fitted 
weighting  function  for  each  channel  and  the  limitation  of  the  functional  form 
used  are  the  two  major  sources  of  errors.  The  latter  problem  is  also  due  to  Che 
limitation  of  the  OMM,  which  is  developed  under  the  assumption  of  an  infinite- 
atmosphere.  To  ameliorate  the  preceding  problems,  measured  radiances  must  btf 
adjusted.  We  define  a  set  of  adjusted  radiances  as  follows: 


Rij  =  ^ijRi 


(3.16) 


where  R,  is  Che  measured  1th  channel  radiance,  and  the  scaling  factor  for  the 
ith  channel  and  the  jth  sharpness  index.  The  scaling  factors  are  introduced  to 
modify  the  measured  radiances  by  Caking  into  account  both  the  effect  of  the 
surface  contribution  and  the  variation  in  the  sharpness  index.  The  proper 
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introduction  of  scaling  factors  are  critical  to  the  success  of  the  retrieval  and 
can  be  approximately  calcul.^ted  by: 

^ij  =  RiVRi‘  . 

where  denotes  tlie  synthetic  radiance  computed  from  Eq .  (3.2),  in  whicJi  the 
coefficients  are  those  determined  in  the  forward  analysis  for  the  ith  channel  and 
the  jth  sharpness  index,  and  R.  is  the  synthetic  radiance  computed  from  Eq . 
(3.13),  using  the  Planck  intensity  profile  for  the  ith  channel  and  weighrinj' 
function.  In  this  way,  the  irregular  behavior  of  the  PH  function  for  p  >  Ps  and 
the  variation  in  the  sharpness  index  can  be  removed  in  the  retrieval  analysis. 
Note  that  also  depends  on  the  spectral  properties  and  atmospheric  temperature 
profile.  Table  7  lists  the  values  of  based  on  the  spectral  properties  of  the 
HIRS  channels  and  the  US  Standard  Temperature  Profile.  The  fitting  of  the  Planck 
intensity  for  k  -  0.49  (j-1)  is  poor.  Thus  only  the  last  six  sharpne.ss  indices 
were  used  in  the  analysis.  We  notice  that  4>ii  <  1.00  for  i  -  1-3.  This  is 
because  the  PH  function  is  unable  to  simulate  the  stratospheric  Planck  intensity 
profile  for  p  <  pj .  For  i  -  4  and  5.  4>ii  values  are  close  to  1  in  view  of  the 
fact  that  the  stratospheric  and  surface  contributions  to  the  radiances  are  small. 
For  i  -  6  and  7,  surface  effects  are  significant  and  become  dominating  for  j  -= 
2  and  3.  Some  of  the  adjusted  radiances  may  not  be  realistic.  Nevertheless, 
they  can  be  used  in  a  consistent  manner  in  the  context  of  the  OMM  to  obtain 
reasonable  Planck  profiles  at  the  desired  pressure  levels. 

For  the  jth  sharpness  index,  the  adjusted  radiances  Rij  are  fitted  to  the  FH 
function,  using  the  prescribed  values  of  mj  and  m2  .  The  coefficients  rj,  r2 ,  hi, 
and  L2  are  subsequently  computed.  These  coefficients  are  used  to  compute  the 
Planck  intensity  profiles  according  to  Eq.  (3.9).  From  each  Planck  intensity 
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Table  7.  Values  of  the  scaling  factor,  based  on  the  channel  properties 

of  seven  HIRS  channels  and  US  Standard  Temperature  Profile.  The 
channel  and  sharpness  index  are  denoted  by  i  and  j ,  respectively. 


i/j 

2 

3 

4 

5 

6 

7 

1 

0.787 

0.788 

0.818 

0.810 

0.821 

0.824 

2 

0.928 

0.929 

0.947 

0.922 

0.920 

0.915 

3 

0.989 

0.847 

0.938 

0.918 

0.929 

0.936 

4 

1.060 

1.080 

0.989 

1.000 

0.925 

0.959 

5 

1.050 

0.999 

1.080 

1.080 

1.010 

1.050 

6 

0.589 

0.192 

0.931 

0.971 

1.070 

1.040 

7 

0 

0 

0.585 

0.639 

0.964 

0.857 

profile,  a  temperature  profile  can  be  obtained.  The  weighted  mean  temperature 
profile  is  then  computed  from  Eq.  (3.10).  Figure  4  shows  the  retrieved 
temperature  profile  using  the  values  of  listed  in  Table  7  and  the  synthetic 
radiances,  Rj,  listed  in  Table  6.  Comparing  Figs.  4  and  3,  we  notice  that  the 
profile  derived  from  the  adjusted  radiances  is  much  more  improved  in  accuracy 
than  that  derived  from  the  original  synthetic  radiances.  The  rms  error  for  p  > 
34.7  mb  is  1.92  K,  and  for  p  >  250  mb,  it  is  about  1.71  K.  The  maximum  error 
occurs  near  the  tropopause  (-  4.5  K)  due  to  the  failure  of  PH  function  to  fit 
exactly  the  simulated  radiance  minima.  The  preceding  temperature  retrieval 
results  using  the  OMM  appear  to  be  encouraging. 

In  order  to  test  the  method  on  a  variety  of  atmospheric  conditions,  we  also 
performed  retrievals  using  synthetic  radiances  using  tropical  and  sub-arctic 
winter  climatological  temperature  profiles  (McClatchey,  1971).  In  each  case,  two 
sets  of  scaling  factors  are  used.  One  is  the  set  listed  in  Table  7  (standard 
scaling  factor) ,  while  the  other  is  the  set  from  synthetic  studies  based  on 
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Retrieved  temperature  profiles  using  the  adjusted  synthetic 
radiances.  The  scaling  factor  set  is  derived  using  the  US  Standard 
Atmosphere . 


respective  climatological  profiles.  Figures  5a  and  b  show  retrieval  results  for 
tropical  and  sub-arctic  winter  atmospheres,  respectively.  The  input 
climatological  profiles  are  also  shown  for  comparison  purposes.  As  expected, 
results  using  the  climatological  scaling  factors  are  more  accurate  than  those 
using  standard  scaling  factors.  The  rms  error  for  tropical  atmosphere  is  2.36 
K  using  tropical  scaling  factors  for  p  >  250  mb.  It  is  4.77  K  using  standard 
scaling  factors.  For  sub-arctic  winter  atmosphere  we  find  the  rms  errors  of  1.87 
K  and  2.58  K,  using  climatological  and  standard  scaling  factors,  respectively. 
In  the  case  of  the  tropical  atmosphere,  serious  errors  occur  near  the  tropopause 
where  there  is  a  strong  temperature  inversion.  As  mentioned  earlier,  the  PH 
function  in  the  present  form  cannot  properly  simulate  the  minima  in  radiances. 
However,  the  retrieved  temperature  profile  in  the  lower  troposphere  using  the 
scaling  factors  is  quite  reasonable. 

Future  efforts  should  be  focused  on  the  investigation  of  the  physical  basis 
of  scaling  factors,  and  the  development  of  correlations  between  radiances  and 
scaling  factors  before  this  novel  method  can  be  applied  to  real  data. 


38 


Fig.  5(a).  Same  as  Fig.  U,  except  for  the  tropical  atmosphere.  The  scaling 
factor  sets  are  derived  based  on  the  US  Standard  Atmosphere 
(dashod-dot  curve)  and  climatological  profiles  (dashed  curve) . 


Fig.  5(b).  Same  as  Fig.  5(a),  except  for  the  sub-arctic  atmosphere. 
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It  follows  that  B^(p)  -  0,  for  p  >  p,.  This  indicates  that,  in  order  to  fit 
B^(p) ,  we  must  find  a  functional  form  that  approaches  zero  for  p  >  Pj.  Many 
mathematical  functions  are  qualified  for  this  purpose.  However,  for  the  purpose 
of  minimizing  error  propagation,  we  shall  use  the  weighting  function  as  our 
functional  form  (Rodgers,  1976).  The  reduced  Planck  intensity  may  be  expressed 
in  terms  of  the  weighting  function  in  the  form 

N 

V  (P)  =  j^i  bj  W5(icj,pj,p)  (3.19) 

The  advantage  of  using  the  weighting  function  as  the  base  function  for  the 
Planck  profile  is  that  the  propagation  of  errors  in  the  radiances  into  the 
retrieved  Planck  profile  is  minimized.  This  is  demonstrated  below.  If  we 
express  B„  as  a  linear  combination  of  arbitrary  functional  forms  such  that 

N 

Bj<P)  =  j^i  bj  K(pj,  p)  ,  (3.20) 

for  discrete  values  p^,  from  Eq.  (2.3)  we  have 

N 

R<Pi>  =  j^l  ^ij  ’  (3.21) 

where 

Aij  =  K(pj,  p)  W-(Pi,  p)  dp  .  (3.22) 

However,  since  bj  is  related  to  B-  via  Eq.  (3.20),  we  may  express  B-  in  terms  of 
R.  Combining  Eqs.  (3.20)  and  (3.21),  we  obtain 

B?(P)  =  K(?j,  p)R(^)  ,  23) 
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where  is  the  element  of  the  inverse  of  matrix,  [Ajj],  By  denoting 


Di(p)  .  S  lAij)-iK(pj,  p) 


(3.24) 


we  have 


B-(p)  -  £  Di(p)R(p  )  . 

i*l 


(3.25) 


Thus,  if  we  minimize  Dj^(p) ,  the  propagation  of  noise  in  R(p)  into  B-(p)  would  be 
reduced  to  a  minimum.  There  is  one  constraint,  however.  For  the  retrieved 
Planck  profile  to  be  exact  at  discrete  points,  we  must  have 


£*  Di(p)  W;(pj,  p)  dp  =  fij 


(3.26) 


where  6^^  is  the  Kronecker  delta  function.  By  jointly  minimizing  every  element 
2 

of  Dj^(p),  subject  to  the  constraint  expressed  by  Eq.  (3.26),  the  best  value  of 
Di(p)  may  be  found  using  the  calculus  of  variation  in  the  form 


Di(p)  =  £  W(pj,  p)  , 


(3.27) 


where 


Wi'j  *  W?(Pj.  P)  Wj;(Pi.  P)  dp  , 


(3.28) 


and  [W. .]'^  is  the  element  of  the  inverse  of  matrix,  W, . . 


By  comparing  Eqs.  (3.24)  and  (3.27)  along  with  Eqs.  (3.22)  and  (3.28),  it  is 


clear  that  if 


K(Pj,  P)  *  Wp(pj,  p)  . 


(3.29) 


the  propagation  of  errors  in  the  radiances  into  the  retrieved  profile  will  be 
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minimized,  since  Eq.  (3.29)  satisfies  the  requirement  of  minimization  as  shown 
in  Eq.  (3.27). 
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COMPARISONS  BETWEEN  DIM  AND  OMM 


In  this  section,  we  shall  first  demonstrate  that  the  DIM  and  OMM  for 
temperature  retrievals  are  equivalent.  Preliminary  results  on  the  application 
of  the  DIM  to  real-time  radiance  data  set  obtained  from  polar-orbiting  satellites 
are  then  presented. 

4. 1  On  the  Equivalence  Between  DIM  and  OMM 

The  basic  theory  of  the  OMM  has  been  described  in  Section  2.  The  radiances 
can  be  expressed  in  terms  of  a  polynomial  function  of  p  as  shown  in  Eq.  (2.1). 
The  Planck  function  can  then  be  computed  from  Eq.  (2.7)  using  this  function.  We 
shall  begin  with  Eq.  (2.1)  and  derive  Eq.  (2.7)  using  the  DIM.  Since  the  DIM  is 
basically  a  method  of  linear  transformation,  we  single  out  an  arbitrary  term  of 
the  nth  power  of  p  in  the  radiance  expression  for  the  sake  of  clarity.  We  will 
attempt  to  show  that  the  Planck  function  can  be  expressed  in  terms  of  the  nth 
power  of  p.  We  may  write 


R(p)  =  p" 


(4.1) 


Using  the  OMM,  the  corresponding  Planck  function  is 


Define 


B(p)  =  p"'^  W-(p)  dp 


p  =  e"’'  or  y  =  -fn  p  , 
so  that  Eq .  (4.1)  can  be  rewritten  as 


(4.2) 


(4.3) 
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R(y)  »  e-"y 


According  to  the  theory  of  the  DIM,  the  Planck  function  can  be  expressed  as  the 
linear  sum  of  higher-order  derivatives  of  R  with  respect  to  the  logarithm  of 
pressure  as  follows: 


B(x) 


Z  A. 
j=0 


dJR(x) 

dx'^ 


(^.5) 


where 

X  =  -inp  , 

and  Aj  are  coefficients  obtained  by  the  Laplace  transform.  Since  R  is  an 
exponential  function,  its  higher-order  derivatives  are  given  by: 

^  ,  (^.7) 

dxJ 

On  substituting  Eq.  (A. 7)  into  (4.5),  we  obtain 


where 


B(x)  =  c  e‘"* 


(4.8) 


c  »  S  A,(-n) 
j=0 


(4.9) 


In  order  to  have  an  analytic  expression  for  Aj,  Eq.  (4.2)  may  be  written  as  a 
bilateral  Laplace  transform  as  follows: 


w(-s)  =  e*’' W(y)dy  ,  (4.10) 

where  s  is  a  transform  parameter.  We  then  expand  l/w(-s)  into  a  McLaurin  series 
in  the  form 
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(w(-s)  ]*^ 


(4.11) 


k=0 

with 

Ak  *  [w(-s)]i^i  /  k! 

Comparing  Eq.  (4.11)  with  Eq.  (4.9),  we  find 

c  =  [w(-n)]-i  . 

Using  Eq.  (4.10),  we  obtain 


(4.12) 


(4.13) 


c 


W-(p)  dp 


(4.14) 


Since  e'"*  -  p",  Eq.  (4.8)  is  identical  to  Eq.  (4.2).  We  have  showed  that  the  DIM 
is  equivalent  to  the  OMM.  Further,  we  may  also  compute  the  coefficients  Aj  from 
known  values  of  c  by  a  matrix  inversion. 

In  an  attempt  to  provide  mathematically  more  rigorous  proof  of  the  uniqueness 
of  the  solution  of  either  the  DIM  or  the  OMM,  Liao  (1989)  showed  that  the 
convolution  transforms  arising  from  geophysics,  such  as  Eq.  (2.3),  are  of  the 
Laguerre- Polya  class.  The  special  characteristics  of  the  Laguerre- Polya  class 
is  chat  its  inversion  function  is  the  convergence  limit  of  any  sequence  of 
polynomials  having  real  roots  only.  It  follows  chat  a  unique  solution  can  be 
obtained  for  the  convolution  transform  whose  kernel  can  be  shown  to  belong  to 
such  a  class.  Details  are  given  in  the  form  of  a  short  research  note  in  the 
Appendix . 

In  our  previous  work  (Liou  and  Ou,  1989),  we  used  a  polynomial  in  inp 
(instead  of  p)  to  fit  radiances.  The  retrieved  temperature  profile  based  on  a 
polynomial  in  fnp  (Fig.  7  in  Liou  and  Ou,  1989)  is  more  satisfactory  than  that 
based  on  a  polynomial  in  p  (Fig.  3)  for  p  >  250  mb.  However,  for  p  <  250  mb,  the 
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DIM  also  suffered  large  errors.  Apparently  the  DIM  depends  less  on  the  surface 
contribution.  To  understand  the  reason  as  to  why  the  DIM  works  better  in  the 
temperature  retrieval,  we  perform  analysis  on  the  behavior  of  the  fitted  Planck 
function  for  p  >  Pj. 

By  fitting  the  radiances  using  the  nth  power  in  -in  p,  we  have 

R(-in  p)  =  (-in  p)"  .  (^15) 

The  Planck  function  is  then  given  by 

B(-in  p)  =  Z  fj  (-in  p)"  . 

J*0  (4.16) 

where  fj's  are  coefficients  to  be  determined.  The  Planck  function  is  also  a 
polynomial  in  in  p.  The  Planck  function  profile  of  HIRS  Channel  7  is  fitted  to 
4th-degree  polynomials  both  in  p  and  in  p,  as  shown  in  Fig.  6.  The  fitting  in 
both  cases  are  quite  satisfactory  above  the  surface.  However,  below  the  surface 
the  fitting  in  in  p  is  much  closer  to  an  isothermal  profile  than  that  in  p.  Also 
shown  in  the  figure  is  the  generalized  weighting  function  of  Channel  7,  which  has 
significant  contribution  for  p  >  1000  mb.  For  this  reason,  the  behavior  of  the 
fitting  for  p  >  p,  will  greatly  affect  the  radiance  of  Channel  7.  The  preceding 
analysis  explains  why  the  DIM  performs  better  than  the  OMM  near  the  surface. 

A . 2  Preliminary  Applications  of  the  DIM  to  HIRS  Data 

We  have  shown  that  the  OMM  needs  modifications  in  order  to  be  practical  in 
tv:mperature  retrievals  and  that  the  DIM  is  better  than  the  OMM.  We  have  made  an 
attempt  to  investigate  the  applicability  of  the  DIM  to  the  inference  of 
atmospheric  temperatures  using  real-time  data. 

An  archive  of  3A73  collocated  clear-sky  temperature  profiles  and  radiance 
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Fig.  6.  The  fitting  of  4th-order  polynomials  in  pressure  (dashed  curve)  and 
in  logarithm  of  pressure  (dashed-dot  curve)  to  the  Channel  7  Planck 
Intensity  profile  (solid  curve)  using  the  US  Standard  Atmosphere. 
Also  shown  Is  the  Channel  7  weighting  function. 
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data  sets  (kindly  provided  to  us  by  Dr.  L.  McMillin)  were  used  to  perform  the 
retrieval  analysis.  The  temperature  profiles  were  based  on  local  sounding,  while 
the  radiance  data  were  obtained  from  HIRS  instruments  on  board  NOAA  10  polar- 
orbiting  satellites  .  For  each  prof lie ,  the  latitude,  longitude,  time,  date,  land 
or  sea,  solar  zenith  angle,  and  water  vapor  mixing  ratio  profile  are  also 
specified.  The  data  are  randomly  distributed  over  lands  and  oceans.  There  were 
1964  cases  over  lands  and  1509  cases  over  oceans.  The  data  were  largely  obtained 
between  March  21  and  April  10,  1987.  All  the  cases  correspond  predominantly  to 
clear  atmospheres  with  a  minimum  amount  of  cloud  contamination.  Temperature 
soundings  over  land  were  obtained  from  land  stations.  Sounding  over  the  oceans 
wore  obtained  from  coastal  stations,  island  stations,  platforms,  and  ships.  The 
observed  radiances  are  corrected  to  the  nadir  direction. 

To  perform  retrieval  exercises,  we  first  convert  the  brightness  temperature 
values  of  Channels  1-7  to  radiances,  we  then  follow  the  procedures  outlined  by 
Liou  and  Ou  (1989)  to  obtain  the  retrieved  temperature  profile.  The  retrieved 
results  are  compared  with  co-located  sounding  profiles.  Specifically, 
temperature  values  at  the  peak  pressure  level  of  the  weighting  functions  are 
selected  and  analyzed.  In  general,  the  rms  errors  at  each  pressure  level  is 
about  the  same  magnitude  as  the  error  from  the  previous  synthetic  computations 
(Liou  and  Ou,  1988).  The  rms  errors  at  the  lower  four  pressure  levels  are  on  the 
order  of  2-5  K.  The  rms  errors  at  the  lowest  pressure  level  (1000  mb)  are  about 
J  K.  The  cause  for  this  error  is  produced  by  the  problem  of  the  surface 
discontinuity . 

Figure  7  shows  the  retrieved  temperature  vs.  sounding  temperatures  at  1000 
mb.  There  are  about  728  points,  which  were  all  located  between  30°N  and  60°N. 
The  straight  line  denotes  that  the  two  temperatures  are  equal.  Most  of  the  data 
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Fig.  7.  The  retrieved  temperatures  vs.  the  sounding  temperatures  at  1000 
mb.  The  radiance  data  were  collected  from  NOAA  9  during  the  period 
between  March  21  to  April  10,  1987,  and  distributed  in  the  latitude 
zone  between  30°N  and  60°N.  The  straight  line  represents  that  the 
retrieved  temperature  is  equal  to  the  sounding  temperature. 
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points  fall  close  to  the  line.  The  correlation  coefficient  is  0.93.  The  errors 
of  its  individual  points  can  be  attributed  to  the  combined  effects  of  many 
factors.  These  include  the  topographic  effects,  cloud  contamination,  uncertainty 
of  the  water  vapor  contamination  in  the  radiances  of  Channel  6  and  7,  errors  in 
the  weighting  function  due  to  different  temperature  profiles,  errors  generated 
due  to  a  least-square  fitting  of  weighting  functions,  errors  in  the  fitting  of 
radiances,  etc. 

In  order  to  improve  the  performance  of  the  DIM,  it  is  necessary  to  undertake 
modifications  and  refinements  of  the  methodology  similar  to  chose  described  in 
Section  3. 
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Section  ,5 

CONCLUSIONS 

In  this  report,  we  first  present  the  basic  theory  involving  the  OMM.  It  is 
shown  that  the  Planck  function  profile  can  be  directly  related  to  the  upwelling 
radiance  through  polynomial  expansions  In  the  pressure  coordinates.  In  order  to 
perform  temperature  retrievals,  an  optical  number  system  has  been  constructed 
based  on  the  generalized  weighting  function.  A  specific  form  for  the  OMM  has 
been  developed  based  on  the  fitting  of  non-linear  hyperbolic  functions  to  the 
measured  radiances.  The  inferred  Planck  profile  Is  shown  to  be  the  combination 
of  the  optical  exponential  function. 

We  apply  the  OMM  to  the  retrieval  of  temperature  profiles  using  the  HIRS-2 
channel  in  the  15  pm  CO2  band.  The  radiances  are  calculated  based  on  the 
spectral  properties  of  seven  channels  using  the  US  Standard  Atmosphere  profile. 
They  are  subsequently  fitted  to  both  a  5th-degree  polynomial  and  a  combined 
polynomial-hyperbolic  (PH)  function.  Retrieved  temperature  profiles  in  both 
cases  deviate  significantly  from  the  standard  temperature  profile.  The  profile 
derived  from  the  PH  function  is  better  than  that  from  the  polynomial.  A  further 
examination  reveals  that  the  retrieved  temperatures  using  the  PH  function  arc 
colder  by  5-12  K  and  warmer  by  0-5  K  below  and  above  250  mb,  respectively. 

To  investigate  the  reasons  for  large  deviations  in  the  retrieved  temperature 
profile,  we  perform  a  forward  analysis.  A  combination  of  a  linear  function  and 
two  terms  of  the  optical  exponential  function  are  fitted  to  the  Planck  profile. 
For  smaller  sharpness  indices,  i.e.,  broader  weighting  function,  the  fitting  is 
poor  and  produces  large  rms  errors  in  the  retrieved  temperature  profiles.  For 
sharpness  indices  larger  than  2.0,  the  rms  errors  are  on  the  order  of  1. 2-2.0  K. 
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In  the  forward  analysis,  we  also  identify  two  problems,  which  are  the  major 
reasons  for  the  deviation  of  the  retrieved  temperature  profile.  The  first 
problem  is  associated  with  the  variation  of  the  sharpness  index  with  the  sensor 
channels.  Fitting  all  the  radiances  to  the  PH  function  using  only  one  set  of 
coefficients  is  mathematically  Inconsistent.  This  problem  can  be  seen  from  the 
variation  pattern  in  the  coefficients  of  the  PH  function  for  different  sharpness 
indices.  The  second  problem  is  caused  by  the  Irregular  behavior  of  the  fitted 
Planck  profile  for  p  >  Ps,  using  the  optical  exponential  function.  The  retrieved 
temperatures  suffer  large  errors  because  the  below-surface  contribution  could  not 
be  properly  accounted  for.  The  below-surface  contribution  is  particularly 
significant  for  channels  with  the  weighting  function  peaks  in  the  lower 
atmosphere.  In  view  of  the  above,  certain  modifications  on  the  OMM  or 
adjustments  to  the  measured  radiances  are  necessary  in  order  for  this  method  to 
be  a  practical  scheme  for  temperature  retrievals. 

We  have  developed  a  scheme  to  adjust  and  augment  the  measured  radiances 
through  a  set  of  scaling  factors,  which  includes  the  effects  of  the  sharpness 
index  variation,  surface  discontinuity,  channel  characteristics,  and  functional 
forms.  We  show  that  mathematically  consistent  retrievals  can  be  performed  on  the 
adjusted  radiances  to  obtain  an  optimum  functional  form,  which  is  most  effective 
in  simulating  the  Planck  profile  and  minimizing  the  error  propagation. 

Through  a  simple  derivation  following  the  fundamental  principles  of  the  DIM, 
we  show  that  the  DIM  and  OMM  are  mathematically  equivalent.  A  rigorous  proof  on 
the  uniqueness  of  the  solution  of  both  methods  is  also  given.  However,  in  our 
previous  work  (Liou  and  Ou,  1988),  we  demonstrated  that  retrieved  temperatures 
based  on  synthetic  computations  using  the  DIM  are  within  about  3  K.  However,  in 
the  present  work,  the  retrieved  temperature  results  using  the  OMM  without 
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adjustments  are  less  satisfactory.  The  superiority  of  the  DIM  is  due  to  the  fact 
that  it  was  the  radiance  polynomial  fitting  in  in  p,  whereas  the  OMM  uses  that 
in  p.  The  fitting  in  the  in  p  coordinate  accounts  for  some  surface  contributions 
to  the  radiances.  Finally,  we  apply  the  DIM  to  an  archive  of  3473  collocated 
temperature  profiles  and  HIRS  radiance  data  sets  from  NOAA  9.  The  accuracy  of 
the  retrieved  temperature  profiles  is  within  about  3  K,  which  is  comparable  to 
that  obtained  from  synthetic  analyses. 
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APPENDIX 


A  Class  of  Laguerre- Polya  Kamel  and  Application 
to  Atmospheric  Temperature  Retrieval 
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ABSTRACT 

A  class  of  convolution  transforms  (generalized  Laplace  transforms)  arising 
from  geophysics  are  shown  to  be  Laguerre -Polya,  l.e.,  their  Inversion  function 
is  the  limit  of  real  polynomials  whose  roots  are  all  real.  It  follows  that,  by 
applying  Widder's  theory,  the  convolution  transforms  can  be  inverted  by  using 
operational  calculus. 
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Section  A.l  Generelixed  Laplace  Transforms 


The  transforms  that  we  will  study  are  of  the  form 

f(x)  -  N(x  -  y)  g(y)  dy  , 

where  the  kernel 


N(x)  »  e*  •  ,  m  >  0 


(A.l) 


To  link  Eq.  (A.l)  to  the  atmospheric  temperature  retrieval,  we  refer  to  King 
(1985).  The  mathematical  problem  is  to  invert  Eq.  (A.l),  i.e.,  to  determine  ^ 

g(x) ,  given  f (x) . 

The  following  analysis  shows  that  if  m  -  1,  Eq.  (A.l)  reduces  to  the 
unilateral  Laplace  transform,  which  Is  usually  given  by 

^^>(3)  =  e'**’ 7(t)  dt  .  (A. 2) 

To  see  this,  let  s  -  e*,  and  t  -  e'^ .  Then  dt  -  -e*^  dy,  and 

Vi(e*)  J  e'**  •  •  •  7(e‘^)  •  e'^  dy  •  (A.  3) 

Thus 

e*  V(e*)  =  f  e*'y  •  •  7(6'^)  dy  . 
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Section  A. 2  Laguerre-Folya  Class 


f 


Definition:  The  reciprocal  of  the  bilateral  Laplace  transform  of  the  convolution 
kernel.  N.  is  called  the  inversion  function  of  the  convolution  transform. 

To  compute  the  inversion  function,  E(s) ,  of  Eq.  (A.l),  we  take  the  Laplace 
transform  as  follows: 

1  .  r“  a*  o-me*/"*  Hv 

■ETiy  J- 

=  J”  m"® 1  u”-""  -  1  e-"  du  .  (^  4) 

by  variable  transformation  u  -  m  e*^”.  Hence 

—  -ms  -  m  »  1 

“  rirnd  -  sH  • 


In  particular,  if  m  -  1,  then 


Since 


E(s) 


1 

ra  -  s-y 


1 

r(s)  r(l  -T5 


sinws 

irs 


(A.  5) 


« 


and 


r(s)-i  =  e^*  s  n 
k=l 


we  have  the  well-known  product  expansion 
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E(s) 


TTT 


iB/k 


where  y  is  Euler's  constant: 


lim 
'  n-^ 


Li.  An  entire  function.  E(s).  is  said  to  belong  to  the  Laguerre-Polya 
class  if  It  can  be  expressed  as 


E(s)  ■  a  •  exp(bs  -  cs^) 


(A. 6) 


where  a.  b.  c.and  Aw  are  real.  C  >  0.  o  Is  a  non-negative  inteEcr.  and 


«0 

s 

k-1 


<  «o 


Clearly,  E(8)  -  Laguerre-Polya. 

The  Inversion  problem  could  be  solved  if  the  inversion  function,  E(s). 
belongs  to  the  Laguerre-Polya  class  as  follows  (Widder,  1971).  Replace  s  in  E(s) 
by  the  differentiation,  D,  to  get  a  differential  operator,  E(D) .  The  function 
g(x)  can  be  determined  by 

g(x)  -  E(D)  f(x) 

in  the  sense  of  operational  calculus. 
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Section  A. 3  The  Inversion  Function 


e 


f 


As  computed  in  Section  A.  2,  the  inversion  function  E(s)  for  Eq .  (A.l)  is 

given  by 


E(s) 


m 


-m*  ♦  ■  -  1 


m  >  0 


(A. 7) 


We  shall  now  confirm  that  E(s)  belongs  to  the  Laguerre-Polya  class. 

Our  argument  is  based  on  the  criterion  described  in  Sheil-Small  (1989,  Lemma 
4),  which  says  that  a  real  entire  function  f(z)  is  in  the  Laguerre-Polya  class 
if,  and  only  if,  ImL(z)  <  0  for  z  t  H* ,  where  L(z)  -  f ' (z)/f (z)  and  H*  is  the 
upper  half  plane. 

Proposition;  Im  <  0  for  s  f 

Proof : 


m""’®(-m)lnm  ra‘“*r'(ra(l  -  s))(-ra) 
“  nma"~s)T  ■  - r(m(l  -  s))2 - 


^  m-">»(-m) 

rWl  -  s)) 


r'Cmd  -  s))' 

r<m(l  -  s)) 

/ 


Let  s  -  X  +  iy,  y  >  0.  We  have 


(A. 8) 


r^(in(l  -  s))  ,  _  1  ^  I  rod  -  s) 

r(ra(l  -  s))  rod  -  s  )  nxl  n(m(l  -  s)  +  n)  ’  (A. 9) 


where  7  is  the  Euler  number.  The  Infinite  series  on  the  right  side  of  the 
equation  is  absolutely  an  uniformly  convergent  in  any  closed  region  of  H"*^. 
Substituting  Eq .  (A. 9)  into  Eq .  (A. 8),  we  get 
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Compute 


and 


Thus , 


1'/"?  -  -n. 

In  m  - 

» 

1  * 

^  *  y. 

'  • 

m(l  -  s) 

E(s) 

mCl  -  s)  n*l 

n(mCl  -  s)  +  n) 

4  . 

-m 


In  m 


m(l  -  s) 


E 

n*l 


m(l  -  s) 
n(m(l  -  s)  ♦  n) 


(A. 10) 


* 


1  =  ^  -  X  ♦  iy 

1  -  S  1  -  X  -  iy  (1  -  x)2  ♦ 


T«  1  -  s  1  -  X  -  iy 

la  —y, - ^ -  ■  la  - i, — 

m(l  -  s)  ♦  n  m  +  n  -  xm  -  imy 


Im  ~  X  -  iy)(P  ♦  n  -  xm  »  imy) 
(m  +  n  -  xm)^  +  m^y* 

t>  (1  ~  x)  my  -  y(m  n  -  xm) 

(m  ♦  n  -  xm)^  +  m^y^ 


my  -  mxy  -  my  -  yn  mxy 
(m  ♦  n  -  xm)^  ■  m^^ 


_ yn 

(m  ♦  n  -  xm)*  ♦  m*y^ 


Im  .  Im 

6(s) 


E 


1  -  s 


1  -  s  n=l  "  a(l  -  s)  '+  n 


_ y 

(1  -  x)^ 


*  E  “ 


-yn 


P  n»l  **  (m  +  n  -  xm)^  ♦  m^y^ 
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